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We analyze the spectrum of the 3-site Bose-Hubbard model with periodic boundary con- 
ditions using a semiclassical method. The Bohr-Sommerfeld quantization is applied to an 
effective classical Hamiltonian which we derive using resonance normal form theory. The 
derivation takes into account the 1:1 resonance between frequencies of a linearized classical 
system, and brings nonlinear terms into a corresponding normal form. The obtained ex- 
pressions reproduce the exact low-energy spectrum of the system remarkably well even for 
a small number of particles N corresponding to fillings of just two particles per site. Such 
small fillings are often used in current experiments, and it is inspiring to get insight into this 
quantum regime using essentially classical calculations. 



I. INTRODUCTION 



Recent experimental and theoretical progress in the field of ultracold quantum gases has stim- 
ulated many studies at the interface of traditionally different disciplines such as atomic physics, 
quantum optics and condensed matter physics providing an intriguing link to fundamental many- 
body problems A specific, but yet particularly interesting topic is the quantum-to-classical 
correspondence in degenerate quantum gases such as the adiabaticity versus non-adiabaticity of 
the quantum dynamics of finite matter- wave systems |3l4l4i| . 

The Bose-Hubbard model (BHM) which we study here is a hallmark of condensed matter theory, 
and was realized experimentally using ultracold quantum gases in optical lattices W following an 



ingenious theoretical suggestion .15|]. At the same time, in the classical limit it is described by the 
celebrated Discrete Nonlinear Schrodinger Equation (DNLSE), which possesses both a rich statics 
and dynamics [3]. Properties of the BHM at high fillings (many particles per site) are known to be 



well-reproduced by several semiclassical methods, such as the mean-field approximation [17|1 or the 
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truncated Wigner approximation jisl. Il9|. At low fillings, one would generally expect semic 
methods to be inapplicable. Interesting enough, in a two-site BHM it was shown recently 
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that semiclassical quantization reproduces the quantum spectrum remarkably well. That is, when 
considering the classical limit of the BHM and applying e.g. the Bohr-Sommerfeld quantization to 
the classical action, one obtains (semi) classical energies which show a good quantitative agreement 
with the exact (quantum) energies. A similar approach has been taken for the case of the three- 



and five-site BHM [2^, with significant insights into the qualitative properties of the quantum 
spectrum, but no explicit expressions for the semiclassical spectrum having been obtained, and a 
corresponding comparison of the quantum and semiclassical predictions is still missing. 

As we show here, the multi-site BHM exhibits an intricate classical dynamics which renders the 
construction of an effective Hamiltonian, necessary for a corresponding semiclassical quantization, 
a nontrivial mathematical problem. Our approach for the derivation of an effective classical Hamil- 
tonian relies on resonance normal form theory [23|]. We note that the classical system linearized 



around its equilibrium can be represented as a collection of harmonic oscillators, the frequencies of 
these oscillators being doubly degenerate. A pair of oscillators with frequencies in resonance 1:1 to 
each other can be analyzed using normalization techniques. To this end one needs to apply a series 
of canonical transformations that bring the quadratic and quartic terms of the Hamiltonian to a 
normal form, thereby completely eliminating all cubic terms. The obtained normalized Hamilto- 
nian, written in terms of action-angle variables, depends on a pair of classical actions, which can 
then be straightforwardly quantized. 

We note that our use of the normalization technique is very similar to that exploited recently 
in the studies of the mean-field dynamics of a nonlinear Stimulated Raman Adiabatic Passage 



(STIRAP) process [24l . l25l |. In its simplest version, nonlinear STIRAP considers three uniform 
condensates: an atomic Bose-Einstein condensate (BEG) in its ground state and a molecular BEG 
in its ground and an excited state coupled by a laser fields (a theory for the non-uniform case has 



also been developed 



261]). Using a certain time-dependent sequence of laser pulses, it is possible to 



convert the atomic BEG into a ground-state molecular BEG without populating the molecular BEG 
in the excited state. Prom a mathematical point of view, the system is a Hamiltonian dynamical 
system with two degrees of freedom possessing a rich dynamics, including e.g. nonlinear instabilities 
due to 1:1 resonances. It is an interesting fact that exactly the same degeneracy influences the 
dynamics of small Bose-Hubbard chains, as we will show below. 

In detail we proceed as follows. Section [H] describes our classical and quantum model. It 
contains a derivation of the effective classical Hamiltonian, its semiclassical quantization and a 
comparison with the exact numerical results for the spectra for both the weak and strong interaction 
regime. Section IIIII provides our conclusions. The Appendix contains a discussion of the two-site 
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Bose-Hubbard chain, where the difference between weak and strong interaction regimes becomes 
transparent. 

II. CLASSICAL AND QUANTUM MODEL 

Let us consider the 3-site Bose-Hubbard model with periodic boundary conditions i.e. a ring 
geometry 

U ^ 

H = -J'^ialak + h.c.) + — "^alalaiai, (1) 

{i,k) 1=1 

where the first sum runs over all the nearest neighbours. As discussed e.g. in Ref. I22I, one 
may try to understand its spectrum using the quantization of a corresponding effective classical 
Hamiltonian which will be derived in the following. To apply a semiclassical quantization, we firstly 
slightly modify the BHM Hamiltonian, thereby writing it in a symmetrized form: 

h=-jy: (aia, + h.c.) + 1 5:(n, + (2) 

{i,k) 1=1 

where hi = a\ai is the particle number operator for the / — th site. Since the total number of 
particles is constant, this modification introduces only a (uniform) shift of all energy levels. The 
classical limit is then obtained by introducing c — numbers for the operators hi in analogy to a 
classical coherent state formalism hi + ^ — )• |'i/'zP- This leads to the Discrete Nonlinear Schrodinger 
(DNLS) equation 



H = -jY^ii^t^k + C.C.) + m\ (3) 
{i,k) 1=1 

'^ = ^ = -J(^,+i + + U^^.,\^l.,\^ (4) 

It is important to note that the normalization of the classical amplitudes is J2i iV'iP = = 
N + ^, where is the total number of atoms in the quantum model [3]. Introducing real pairs 
{xn,yn) with tpn = , the DNLSE is equivalent to the Hamiltonian equations of motion of the 

Hamiltonian 

3 ( 2 , 2\2 

H = UY. _ J ^ (^^^, + y^y,) , (5) 

1 = 1 {i,k) 
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Switching to polar coordinates Xi = \/2nism 4>i, Ui = \/2nj cos we arrive at the Hamiltonian 

U ^ 

H = —^nf - 2J^ {^/r^n^cos{(j)i - (pk)) , (6) 
1=1 (i,k) 

with J2i ''T'i = const = Ns = N + ^. We now perform a rescaling Ui = liNg, H = hNs. With 
g = NgU this leads to the Hamiltonian 



3 

^= f E^'-2«^E V^cos(0,-0fc), (7) 
1=1 {i,k) 

with h = const = 1. 

It is possible to introduce an effective classical Hamiltonian for the low-energy and high-energy 
dynamics of the system. Here we restrict ourselves to the low-energy part. 

For g > gcr = —9/2 the ground state is uniform in density and phase /j = 1/3 Vi and (pi = 
Vi, J. For strongly attractive interaction g < gcr = —9/2 the ground state is essentially different 



16|. Here we only consider values of interactions far from this bifurcation, i.e. we have g > gc 



A. Low-energy spectrum for small and moderate values of interaction 

It is possible to eliminate one degree of freedom by applying a transformation with the generating 
function W = pi{(j)i - </>2) + P2{4>1 + (t>2 - '^4>?.) / + Pz{(pi + 02 + h)/^- 

This generating function is chosen such that the expressions occurring in the below-given Taylor 
expansion will have a simple appearance. 

The transformed Hamiltonian depends only on two new phases 6i , 02 and their corresponding 
momenta pi,P2) while the third momentum is an integral of motion = Y^^=i h = ^- At values 
of g larger than the critical value gcr, the stable equilibrium is at the origin pi 2 = 0, 9i^2 = 0. We 
expand our Hamiltonian around the origin up to terms of fourth-order power in coordinates and 
momenta. Without loss of generality we put J = 1 in the following. 

The resulting Hamiltonian is 

h = ho + ^^^ + '^-±^{pl+pl) + H, + H,, (8) 

where and Hi^ contain cubic and quartic terms, respectively, and Kq = —1 -|- | is a constant. 
Note that the quadratic part of the Hamiltonian describes two harmonic oscillators whose frequen- 
cies are in resonance 1:1. It is exactly this type of degeneracy which was recently considered in 



Ref. 



24 



25l |. The quadratic part of the Hamiltonian implies a bifurcation at gcr = —9/2, i.e. at 
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sufficiently strong attractive interaction. As already mentioned above, we focus on the low-energy 
dynamics off this bifurcation, i.e. for weakly attractive or repulsive interactions. 
The cubic terms are given by 

/3 

H3 = Y + {el - el - 27pl)p2 + 9pl) (9) 

To bring the Hamiltonian to its normal form, we need to get rid of the cubic terms. This can 
be done by a nonlinear near-to-identity canonical transformation pi 2, 61,2 Pi,2j^i,2 determined 
by the generating function 



W = ei (Pi + PP1P2) + ^2 {P2 + aPi + 7Pf + cdl + del), (10) 
with the coefficients 

a = -^^, /0 = -2q;, 7 = -", (H) 
18 — ^ 

"^^^^^ d = -c/3, ^ = 9 + 25 (12) 

This transformation does not change the quadratic terms, removes the cubic terms, and mod- 
ifies the quartic ones. We subsequently change to polar coordinates Ji,2,'5'i,2, with Xi^2 = 
V2Ji^^i/^sin$i,2,Pi,2 = V^Ji^^-^/^cos$i,2 

In the resulting Hamiltonian, we keep only slowly varying terms which depend on ^1 — ^2, and 
average out (i.e. omit) other ('fast') trigonometric terms, e.g. cos 2$i, cos 2$2, cos(2($i + ^2)), 
etc. We thus arrive to the normal form 



H = Va{Ji + J2) + B{jf + J|) + CJ1J2 + DJ1J2 cos[2($i - $2)], (13) 

where 

B = -^{^^ + ^^9 + 9% C = -^{6 + g){18 + g), (14) 

^ = -i|2(-54 + 245 + <7'), (15) 

One can easily see that at 5 = all coefficients B,C,D are equal to zero, implying the degeneracy 

of the Hamiltonian with respect to Ji — J2 at this point. Introducing J = "^^^"^^ , P = we 
get the Hamiltonian 
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h = 2VAJ + 2B{J^ + P^) + CiJ^ - + D{J^ - cos 2$ (16) 

where P and ^ are canonically conjugate. This Hamiltonian is integrable, since J =const. 
J is the first action variable of our effective Hamiltonian, the second one we find by calculating 
K = ^ J Pd^. This integral can be calculated analytically, which provides us with an expression 
for K{h). Inverting it, we finally find that 



h = 2nj + cJ^ - e[K^ - 2KJ] , (17) 

where 

c = 2C = -^{6 + g){18 + g), e = D = + 2Ag + g^), Q = VA, (18) 

The above expression ()17p together with the constants (jlSp constitute the main result of this 
paper. The quantization of the actions J, K should reproduce the low-lying energy levels of the 
system. We use the following quantization of the actions 

77-1-1 

2J = n = 0,l,..N, (19) 

2K = ^^^±p, m = 0,..n, (20) 

The resulting semiclassical energy-levels are shown in Figs. [1] [2l in comparison with 

the results from exact numerical diagonalization for N=24 and N=6. It is observed that for 
moderate values of the interaction strength \g\ the exact spectrum is reproduced very well. To 
be more specific, there are two phenomena occuring for the spectrum with increasing \g\. At 
zero interaction the exact spectrum coincides with the Bogoliubov one, and the energy levels are 
organized in (degenerate) Bogoliubov bands. As g is increased, these energy bands decrease and 
their degeneracy is lifted. 

The first phenomenon is reproduced by the semiclassical approach remarkably well. That is, we 
see from Figs (TJ [2] that the energetical distance between the exact solutions and the semiclassical 
prediction is much less than the distance between the exact and the linearized solutions (i.e., the 
deviation of the exact spectrum from the Bogoliubov frequencies is much larger than its deviation 
from the semiclassical prediction). At the same time, the spreading of the energy levels within the 
Bogoliubov band is reproduced not very well. For large (^-values the deviation is considerable. As 
g is increased, fiuctuations with respect to the phase grow and an expansion around the classical 
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FIG. 1: Energy levels of the three-site BHM (red solid lines) and semiclassical levels given by Eq. ([T7)) (black 
dashed lines) as a function of the interaction strength g ioi N — 24 atoms. The energy levels are rescaled 
to the Bogoliubov frequencies, i.e. we show {E — Eo)/il, where Eq is the energy of the ground state, and 
fl = \J~A = v^9~F25. The Bogoliubov levels are denoted by blue dotted lines. 

ground state becomes inapplicable: the phase is not confined to the vicinity of anymore. The 
case of large g is analyzed in the next section. 

B. LoviT-energy spectrum for strong interactions 

At strong interactions, even though the classical ground state remains unchanged (uniform in 
density and phase), the semiclassical ground state becomes qualitatively different. The semiclassi- 
cal ground state includes zero-point oscillations, i.e. the corresponding classical trajectories possess 
certain non-zero classical actions. With increasing interaction strength, the area of the phase space 
filled with trajectories oscihating around the classical ground state shrinks. At a certain value of g 
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FIG. 2: Energy levels of the three-site BHM (red solid lines) and semiclassical levels given by Eq. ([TT]) 
(black dashed line) as a function of the interaction strength g ior N = 6. The energy levels are rescaled to 
the Bogoliubov frequencies. The corresponding Bogoliubov levels are denoted by blue dotted lines. 

the oscillatory area cannot accommodate trajectories with these minimal actions. As a result, the 
nature of the classical motion changes from oscillatory to rotational. In other words, the zero-point 
oscillations destroy the phase coherence. This is most easily illustrated for the double-well case 
(see supplementary information). One can consider it as a (semi)classical counterpart of the Mott- 
insulator transition. To be more precise, in finite quantum chains the Mott-insulator transition 
becomes a smooth crossover, and the change of the classical motion described above is the classical 
counterpart of this crossover. 

From a technical point of view, this qualitative change in classical trajectories corresponding 
to the ground and low-lying states leads to the inapplicability of the fourth-order expansion of 
the classical Hamiltonian around the origin: phases now cover the interval [— vr, vr) and are not 
restricted to the vicinity of zero. An adequate effective Hamiltonian will consequently change 



severely. To derive it, one can use classical perturbation theory such as Linstedt's method 
Similar to the case of weak interactions, we introduce 

3 

2 6 



Pi = ^-^^, P2 = \{h+h-2h), P3 = ^/, = l, (21) 

i=l 

Expanding the Hamiltonian in powers of Pi,P2 and 1/g, one obtains 

H = 9[l + P? + 3i^| - ^(cos $1 + 2 cos ^ cos ^) + 0{Pjg)] (22) 

We assume that the interaction strength g is large enough such that Pi > \/l/g holds. This 
holds if 1 /N ^ i/l/g. Then, the classical dynamics takes place outside the resonance zones i.e. it 
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becomes a rotational motion. Then, in the first approximation of Linstedt's method, one introduces 
new actions Pi — t- Ji , P2 — ^ J2-, and the averaged Hamiltonian in the new actions Ji , J2 coincides 
with the unperturbed initial Hamiltonian (this result can be obtained by other methods as well 
0), 

H = g[]: + Jl + 3J|] (23) 
b 

The semiclassical energy levels are obtained by quantizing Ji , J2 : 
„ AT /I 2 2\ ni — n2 ni + n2 1 

Eni,n2 = gNs{-+ jl + j2),jl = , -^2 = ^7^ 3' "'1=0'"'^' "2 = 0, .., iV - m, 

(24) 

FigE] shows a comparison of the semiclassical and exact results. In the second approximation 
of Lindstedt's method, the dependence of the new Hamiltonian on the actions is slightly modified 
and the corresponding expressions are more cumbersome and therefore not provided here. 

III. CONCLUDING REMARKS 

We have shown that by performing a proper analysis of the classical counterpart of a few-site 
Bose-Hubbard model one can gain valuable insights into the dynamics of the system. 

The classical system possesses two degrees of freedom and its small-amplitude oscillations around 
the ground (equilibrium) state can be brought into the form of two decoupled linear oscillators. 
An important feature of the system is that the oscillators are in resonance 1:1 to each other. 
As a subsequent approximation, considering excitations on top of the ground state with larger 
amplitude, one should therefore apply resonance normal form theory. This allows us to bring the 
system to the form of an integrable nonlinear Hamiltonian depending on two classical actions. The 
quantum spectrum can then be reproduced by quantizing these actions. From the point of view 
of the physics of Bose-Einstein condensates, this procedure can be seen as an extension of the 
Bogoliubov transformations to the realm of nonlinear oscillations. 

For strong interactions, the properties of the system change drastically. Here perturbation 
theory with respect to the inverse interaction strength leads to a system of uncoupled rotators. 

An interesting dynamics occurs if one dynamically sweeps the interaction from strong to weak 
values (or vice versa), passing through the crossover region. Depending on the sweeping rate of the 
interaction, the crossover region is passed adiabatically or non-adiabatically and correspondingly 
different amounts of excitations are produced at the end of the sweep. This problem for the triple- 
well is left for future research. However, we solved it for the double-well. This allows, e.g, to 
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FIG. 3: Energy levels of the three-site BHM (black squares) and semiclassical energy levels given by Eq. 
([M)) (red circles) as a function of the interaction strength g. The levels are almost indistinguishable on the 
scale of (a), (b) represents an enlarged view of a part of (a). 

construct a mathematical theory of slow decoupling of two superfluids, extending the results of Q] 
on abrupt decoupling of superfluids to the case of slow sweeps. This question is presented elsewhere 
28|. 

The results can be generalized to the case of longer chains. While it is understood that normal 
form theory is useful for weakly interacting wave dynamics j^l, we are not aware about corre- 
sponding explicit calculations in BHM. 
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Appendix A: Quantization of the two-site chain 

For comparison, we consider here the quantization in case of the double-well problem. The 
two-site Bose-Hubbard Hamiltonian reads 



H = -J{a\ai + a\a2) + ^[ni{ni - 1) + n2(n2 - 1)] (Al) 

The symmetrized form of the interaction term is 

+ \f + (n2 + \?-\- 2(ni + n2)], (A2) 

therefore we will consider below the slightly modified form of the BHM, additionally asuming 
J= 1: 



H = -{alai + a\a2) + ^[{m + + (^2 + ^f] (A3) 
The semiclassical limit is given by the Hamiltonian 

H = -{r2i'i + rM + lov-ii' + m% (A4) 

with the corresponding equations of motion 

dt dij*' dt dijj' ^ ^ 
The normalization of the semiclassical variables is j^] 

|V'iP + lV'2|^ = iV + l = A^s (A6) 
Introducing = , the equations of motion for Xi , yi are determined by the classical Hamilto- 



nian 



H = ?[(^? + vlf + (4 + y^f] - ixiX2 + yiy2), (A7) 
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which after transforming Xi = ■y/2Ii cos (f)i, yi = -i/2Ti sin (f)i becomes 

H = ^{If + I|) - 2y7i72 cos ((/-i - (^2) (A8) 

with the integral of motion I1+I2 = Ng. Employing the generating function W = ^(^1 — ^2) + 
^(01 + ^2) to canonically transform the above Hamiltonian we arrive at 

H = ^{P^ + N^)- - P2 cos 2^ (A9) 

Subsequently the transformation P = Ngp, H = Ngh, and gives us 

h = ^- sjl-p^ cos 24), (AlO) 

where g = UNg. For weak interactions, a phase point oscillates around the ground state. For 
strong interactions, separatrix h = shrinks to the vicinity of p = 0, and when the area within the 
separatrix divided by 27r becomes less than the minimal action 2^ , nature of motion of the phase 
point corresponding to the semiclassical ground state changes: it is a rotation, with (f) covering the 
full interval (0,27r). 

While it is not difficult to calculate the classical action here, we would like to follow a general 

scheme which could be generalized to a M — site Bosc-Hubbard chain. Therefore for weak inter- 
actions we proceed with an expansion of the Hamiltonian around its equilibrium, keeping terms of 
up to the fourth order. Expanding the Hamiltonian close to the origin, one gets 

h = -l + 24>'+p'{l + g/2)/2 + ^ -pV^ - ^ (All) 
The quadratic part can be transformed to action-angle variables easily 



1 



p = y 2//a;cosa;, (p = V2Iuo sin x, io = -yJl + g/2 (A12) 



Then we obtain 



h = -i + ni- + I^Fix), (A13) 
where / F{x) = 0, Q. = Aui = ^2(2 + g). Averaging over x, we get the effective Hamiltonian 

x=0 

h = -i + ni-^-^^i' (A14) 

The quantization of the action leads to /„ = '^'^^^ , n = 0,..N and finally gives us the 
semiclassical energy levels 

E = iV.(?_l + Si,„_?to+M], (A15) 
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FIG. 4: Energy levels of the two-site BHM (black solid lines) and the first 5 semiclassical levels (jAlSp (blue 
dashed lines) for iV = 9 atoms. The first three levels are almost indiscernible on this scale. 

vi^hich reproduce the low-energy spectrum of the dimer BHM remarkably well even for ~ 10 
(see Fig. H]). 
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